# Goal: Make Figure D1
# Dependencies: "reffs_samples.Rdata"

setwd("~/Replication Package")

library(coda)
library(plotrix)

load("reffs_samples.Rdata")

reffs.samples <- as.mcmc.list(reffs.samples)

coef.samples <- reffs.samples[ , substr(names(reffs.samples[[1]][1, ]), 1, 12) != "Sigma.beta[1" & substr(names(reffs.samples[[1]][1, ]), 1, 6) != "beta[1" & substr(names(reffs.samples[[1]][1, ]), 1, 9) != "mu.beta[1"]

coef.samples.all30 <- NULL
for (j in 1:30) {
  coef.samples.all30 <- rbind(coef.samples.all30, coef.samples[[j]])
}

N.iters <- dim(coef.samples.all30)[1]

beta.samples.type2 <- array(coef.samples.all30[ , substr(colnames(coef.samples.all30), 1, 6) == "beta[2"], c(N.iters, 3, 13))
beta.samples.type3 <- array(coef.samples.all30[ , substr(colnames(coef.samples.all30), 1, 6) == "beta[3"], c(N.iters, 3, 13))
beta.samples.type4 <- array(coef.samples.all30[ , substr(colnames(coef.samples.all30), 1, 6) == "beta[4"], c(N.iters, 3, 13))

beta.means.type2 <- apply(beta.samples.type2, c(2, 3), mean)
beta.means.type3 <- apply(beta.samples.type3, c(2, 3), mean)
beta.means.type4 <- apply(beta.samples.type4, c(2, 3), mean)

beta.ql.type2 <- apply(beta.samples.type2, c(2, 3), quantile, p = 0.025)
beta.ql.type3 <- apply(beta.samples.type3, c(2, 3), quantile, p = 0.025)
beta.ql.type4 <- apply(beta.samples.type4, c(2, 3), quantile, p = 0.025)

beta.qu.type2 <- apply(beta.samples.type2, c(2, 3), quantile, p = 0.975)
beta.qu.type3 <- apply(beta.samples.type3, c(2, 3), quantile, p = 0.975)
beta.qu.type4 <- apply(beta.samples.type4, c(2, 3), quantile, p = 0.975)

##-------- FIGURE D1 --------##

jpeg("FigureD1.png", width = 900, height = 600)

par(mfrow = c(1, 3)) 
par(oma = c(4, 10, 0, 0),mar = c(0, 2, 5, 3)) 

varnames <- c("Education", "Age", "Female", "Non-white", "Income proxy", "Ideology (left-right scale)", "National economy", "Personal economy", "Crime victim", "Perception of corruption", "National capital", "Big city")

# (L, H) vs. (L, L)

plotCI(y = 1:12 - 0.1, x = rev(beta.means.type2[1, 2:13]), err = "x", ui = rev(beta.qu.type2[1, 2:13]), li = rev(beta.ql.type2[1, 2:13]), axes = FALSE, ann = FALSE, ylim = c(0.7, 12.1), xlim = c(-1.2, 1.2), pch = 16, pt.bg = par(cex = 0.75), sfrac = 0.005, lwd = 1, col = "black")
plotCI(y = 1:12, x = rev(beta.means.type2[2, 2:13]), err = "x", ui = rev(beta.qu.type2[2, 2:13]), li = rev(beta.ql.type2[2, 2:13]), pch = 16, pt.bg = par(cex = 0.75), sfrac = 0.005, lwd = 1, add = T, col = "gray50")
plotCI(y = 1:12 + 0.1, x = rev(beta.means.type2[3, 2:13]), err = "x", ui = rev(beta.qu.type2[3, 2:13]), li = rev(beta.ql.type2[3, 2:13]), pch = 16, pt.bg = par(cex = 0.75), sfrac = 0.005, lwd = 1, add = T, col = "gray80")
axis(side = 1, cex.axis = 1)
axis(side = 2, las = 1, at = c(1:12), cex.axis = 1, labels = rev(varnames), las = 2)
mtext("Agitator vs. Outsider", at = 0, line = 1, cex = 0.9)
abline(v = 0, col = "grey63")

legend(0.5, 11, c("2014", "2012", "2010"), col = c("gray80", "gray50", "black"), pch = 16, lwd = 1, bty = "n", cex = 1, y.intersp = 2)

# (H, L) vs. (L, L)

plotCI(y = 1:12 - 0.1,x = rev(beta.means.type3[1, 2:13]), err = "x", ui = rev(beta.qu.type3[1, 2:13]), li = rev(beta.ql.type3[1, 2:13]), axes = FALSE, ann = FALSE, ylim = c(0.7, 12.1), xlim = c(-1.2, 1.2), pch = 16, pt.bg = par(cex = 0.75), sfrac = 0.005, lwd = 1, col = "black")
plotCI(y = 1:12, x = rev(beta.means.type3[2, 2:13]), err = "x", ui = rev(beta.qu.type3[2, 2:13]), li = rev(beta.ql.type3[2, 2:13]), pch = 16, pt.bg = par(cex = 0.75), sfrac = 0.005, lwd = 1, add = T, col = "gray50")
plotCI(y = 1:12 + 0.1, x = rev(beta.means.type3[3, 2:13]), err = "x", ui = rev(beta.qu.type3[3, 2:13]), li = rev(beta.ql.type3[3, 2:13]), pch = 16, pt.bg = par(cex = 0.75), sfrac = 0.005, lwd = 1, add = T, col = "gray80")
axis(side = 1, cex.axis = 1)
mtext("Conventional vs. Outsider", at = 0, line = 1, cex = 0.9)
abline(v = 0, col = "grey63")

legend(0.5, 11, c("2014", "2012", "2010"), col = c("gray80", "gray50", "black"), pch = 16, lwd = 1, bty = "n", cex = 1, y.intersp = 2)

# (H, H) vs. (L, L)

plotCI(y = 1:12 - 0.1, x = rev(beta.means.type4[1, 2:13]), err = "x", ui = rev(beta.qu.type4[1, 2:13]), li = rev(beta.ql.type4[1, 2:13]), axes = FALSE, ann = FALSE, ylim = c(0.7, 12.1), xlim = c(-1.2, 1.2), pch = 16, pt.bg = par(cex = 0.75), sfrac = 0.005, lwd = 1, col = "black")
plotCI(y = 1:12, x = rev(beta.means.type4[2, 2:13]), err = "x", ui = rev(beta.qu.type4[2, 2:13]), li = rev(beta.ql.type4[2, 2:13]), pch = 16, pt.bg = par(cex = 0.75), sfrac = 0.005, lwd = 1, add = T, col = "gray50")
plotCI(y = 1:12 + 0.1, x = rev(beta.means.type4[3, 2:13]), err = "x", ui = rev(beta.qu.type4[3, 2:13]), li = rev(beta.ql.type4[3, 2:13]), pch = 16, pt.bg = par(cex = 0.75), sfrac = 0.005, lwd=1, add = T, col = "gray80")
axis(side = 1, cex.axis = 1)
mtext("Activist vs. Outsider", at = 0, line = 1, cex = 0.9)
abline(v = 0, col = "grey63")

legend(0.5, 11, c("2014", "2012", "2010"), col = c("gray80", "gray50", "black"), pch = 16, lwd = 1, bty = "n", cex = 1, y.intersp = 2)

dev.off()
